含教程 | 时间导数分布修复(TDDR)：fNIRS运动校正方法

导读

功能近红外光谱(fNIRS)是一种光学神经成像技术，它作为一种研究皮层活动的工具受到了广泛的关注。由于在头部放置光极，头部运动产生的伪影相对来说比fMRI要小。然而，在数据处理过程中仍需删除运动伪影。研究者提出了一种新的基于鲁棒回归的运动校正方法(即，时间导数分布修复(TDDR))，该方法有效地消除了基线偏移和峰值伪影，而且不需要用户提供任何参数。模拟结果表明，该方法比其他5种运动校正方法具有更好的激活检测性能。在对7-15岁儿童样本的工作记忆任务实证验证中，该方法产生了比其他任何方法更强、更广泛的激活。TDDR校正方法增强了fNIRS作为功能性神经成像方式的可行性，可用于不适应fMRI研究的人群。

前言

fNIRS作为一种研究皮层活动的光学神经成像技术，近年来引起了人们的广泛兴趣。由于fNIRS光电装置放置在参与者的头皮上，fNIRS对头部运动产生的伪影不像功能性磁共振成像(fMRI)等固定传感器方法那么敏感。然而，运动伪影的存在可能仍然是有问题的，特别是在头部运动较多的人群中，如婴幼儿。在测试具有较高头动水平的群体的非典型神经活动模式的研究中尤其如此，如自闭症谱系障碍、注意缺陷/多动障碍和癫痫。如果不加以控制，由头部运动产生的伪影可能在神经活动模式中产生虚假的差异，从而掩盖了真实的差异。因此，尽管fNIRS对头部运动相对不敏感，但仍然有必要采用预处理策略来消除过度头动造成的伪影。

虽然针对fNIRS数据已经开发了几种运动校正算法，但它们都存在一些缺陷。例如，基于小波的方法很好地处理了运动峰值，但加剧了基线偏移伪影。其他方法如运动伪影减少算法(MARA)和目标主成分分析(tPCA)依赖于伪影检测，这需要用户提供几个参数。这些参数在不同仪器、参与者或脑区之间不一定相同，这给用户带来了很大的负担。tPCA和基于相关性的信号改善(CBSI)方法分别对主成分大小、含氧血红蛋白和脱氧血红蛋白之间的关系做出的假设并不总是成立的。考虑到现有方法的这些局限性，本研究的目标是开发一种方法，可以有效地去除这两种类型的运动伪影，无需调优参数，并对被校正的数据进行最少的假设。

在这里，研究者介绍了一种新的伪影校正方法，时间导数分布修复(TDDR)。通过模拟来评估该方法的相对性能，其中任务激活管道在数据上执行，这些数据包括：1)无运动伪影，2)未校正的运动伪影，或3)使用TDDR方法或5种比较方法之一校正的运动伪影。为了在真实数据上证明这种方法的可行性，研究者在儿童样本中进行了一项实证研究，因为该群体通常比成年人群体有更高的运动。研究者让7-15岁的儿童参与者执行了工作记忆任务，因为大量的fMRI和fNIRS研究一致表明，工作记忆能够有效激活儿童前额叶皮层，使其成为验证本研究算法的理想认知过程。然后，使用TDDR方法评估任务激活，并与5种现有的运动校正方法进行比较。

方法

理论

测量的fNIRS信号(x)可以表示为时间(t)的函数，其中信号是所有先前实时波动(y)的时间积分。

这些波动包括反映认知或运动处理的血流动力学反应(主要感兴趣的信号(h)，与系统(生理)调节相关的血流动力学活动(例如，心脏和呼吸振荡和迈耶波)；这些信号有助于记录来自头皮和头骨等组织表层的信号)(p)、运动伪影(m)和仪器噪声(ε)信号。

由于仪器噪声具有随机性和宽频性，在预处理中很难用任何分析方法去除，因此必须在硬件层面上进行处理。本研究假设ε_t可以忽略不计。由于系统生理信号是相对窄带和拟稳态的，因此可以用许多现有的校正方法去除它们，如独立成分分析。此外，由于这些信号主要是来自表层的信号，它们可以通过放置在靠近源的其他探测器记录下来，然后从感兴趣的信号中进行回归。该算法旨在减少运动伪影(m_t)。这里概述的方法依赖于三个假设：

①非运动波动(即非运动相关活动的导数)近似正态分布。

②大部分波动不包含运动伪影。

③当运动伪影存在时，与运动伪影(即它们的导数)相关的波动比非运动波动的幅度大得多。

如果这些假设得到满足，来自运动的信号具有振幅大且不频繁的特点，因此驻留在信号波动的正态分布的遥远尾部。解决该问题的一个可能方法是减少异常大波动的权重。

其中w是范围为[0，1]的权重函数。一个明显的选择是使用基于Chauvenet准则的固定阈值：

当信号波动的幅度大于均值(μ)的α个标准差(σ)时，将其设置为零，其中α通常为3-5左右的某个值。虽然可以通过这种方式为包含异常值的个体波动的分类设置一个界限，但其实也带来了阈值选择的负担，而阈值选择可能因人群、个体或脑区而异。相反，研究者提出了一种类似鲁棒回归中使用的迭代重加权方案。

最常见的鲁棒估计量是m估计量(指所有用通过找到最大化/最小化某个目标函数的参数的总称，也叫极值估值法)。这些估计量被定义为使表达式最小化：

其中d_i是预测值和观测值之间的尺度偏差，ρ是一个权重函数。m估计量是极大似然估计量的概括，最小二乘解是一个特定的非鲁棒情况：。在鲁棒m估计的情况下，权重函数将较低的权重放在较大的偏差上。一个常见的鲁棒估计是Tukey双权函数：

d_t是每个观测值的比例残差。该函数将较低的权重放在远离均值的观测值上，极端偏差减少到零。对于该估计量和许多其他估计量，有必要采用迭代重加权方案，其中参数、残差和观测权重都要迭代重计算，直到参数估计收敛于一个值。使用这种迭代重加权方法，可以产生一组权重，以缩小通常与头部运动相关的过大波动(例如，|d_t|≥1)。

实际考虑

图1A展示了这种方法在真实数据中的效果。由于TDDR利用了信号的时间导数，因此它基本上取决于该信号的采样率。这是因为采样周期与时间导数值的大小直接相关。此外，该算法依赖于对时间导数方差的估计来计算鲁棒权值。由于这个原因，高频(HF)成分的存在(无论是由于仪器/测量噪声还是生理原因)夸大了估计的方差，从而降低了其有效性。为了解决这些问题，本研究首先检查采样率是否足够高，以至于会出现过多的高频成分(＞1Hz)。如果是这样，则应用0.5Hz截止的低通(LP)滤波器并保留滤波器残差，将信号分为低频和高频成分(图1B)。然后将TDDR算法应用于信号的低频成分，再将其添加回未校正的高频成分。从图1C中可以看出，HF成分的时间导数的方差远远大于即使是较大的运动伪影。这种LP滤波过程有效地缓冲了核心算法中采样率过高和高频成分存在的问题。

图1.在真实数据上演示TDDR。

虽然这种算法本质上是迭代的，但它仍然非常快。这是因为每次迭代只涉及到计算信号导数的加权平均值，然后计算新的权重集。所有操作在计算上都相对简单。为了说明这一点，一个模拟的10min扫描，32个通道以20Hz采样，在2.40GHz工作站上以单线程模式应用TDDR大约需要0.13s。在本研究的模拟和实验数据中，权重通常在不到20次迭代中收敛。实现TDDR算法的源代码可以在https://github.com/frankfishburn/TDDR上访问。

模拟验证

数据集生成：模拟的fNIRS信号是使用一个随机生成的自回归模型。32个通道的信号以20 Hz的频率产生，空间协方差为0.33。任务结构包括10min的扫描，再细分为固定时长的block，随机选择1到60s，block之间的间隔为10-60s。通过将任务boxcar与典型血流动力学响应函数进行卷积得到模拟的任务响应。然后将任务响应添加到一半通道的自回归噪声中，振幅为噪声标准差的0.5-25倍。通过产生正弦波，将系统生理振荡添加到由此产生的血红蛋白信号中，正弦波的频率来自心脏(mu 1，sigma 0.1)、呼吸(mu 0.25，sigma 0.025)和迈耶波(mu 0.1，sigma 0.01)噪声的正态分布。使用标准差为0.1 rad/s的正态分布将累积相位漂移添加到每个信号。然后将这些血红蛋白浓度信号转换回原始强度。然后以每分钟2个峰值和2个基线偏移的频率生成峰值和基线偏移运动伪影。峰值被建模为指数增长曲线(x(t)=b^t/τ)，从标准差为25倍输入信号标准差的正态分布中随机选择一个正增长因子(b)，并选择时间常数(τ)产生最终峰值持续时间0.1-10.0s。基线偏移被建模为随机的正或负偏移，振幅选择的方式与峰值相同。以这种方式生成了10000个文件。使用NIRS Brain AnalyzIR工具箱进行模拟。图2显示了一个没有伪影、未校正伪影和通过TDDR校正伪影的模拟时间进程样本。

图2.在模拟数据上演示TDDR。光密度显示为无运动伪影信号(蓝色)，未进行伪影校正的相同信号(红色)，以及使用TDDR校正伪影后的信号(绿色)。可以看到，校正后的信号与无运动伪影信号大致相同。

分析#1：采样率和滤波的影响

预处理。原始强度信号被转换为光密度。然后将信号重采样至120Hz。然后，重采样的光密度信号要么直接应用TDDR，要么使用TDDR只校正信号的低频部分。然后使用修正的Beer-Lambert定律将校正后的光密度信号转换为血红蛋白浓度。为了避免采样率直接改变激活检测的性能，然后将所有信号重采样至1Hz。

激活。使用GLM评估激活程度。为了考虑信号中的低频漂移，模型中包含了离散余弦变换矩阵项，最大频率为1/128Hz。使用受试者工作特征(ROC)曲线评估激活检测的性能。为了在ROC曲线上生成边界，10000次迭代(总共32000个通道)的结果被划分为100个bin，并在曲线上计算均值和标准差。对于全局性能指标，曲线下面积(AUC)通过梯形积分计算。为了验证所选择的迭代次数，每个测试的结果(例如，每个采样率和预滤波)被随机分成100次，AUC中的均方误差(MSE)用计算。每次检验的AUC均方误差均小于10^-5。

分析#2：方法比较

预处理。模拟的fNIRS信号使用TDDR和5种不同的比较方法进行校正，在Homer2中实现。然后根据修正的Beer-Lambert定律将信号转换为血红蛋白浓度，重采样至1Hz。

比较方法。基于相关性的信号改善(CBSI)依赖于以下假设：氧合血红蛋白和脱氧血红蛋白呈强负相关，而运动伪影在这些测量值之间呈正相关。本质上，脱氧信号被缩放到与氧合信号相同的方差，然后从氧合信号中减去。其基本原理是，信号的正相关特征(即运动)将被减法抵消，而负相关特征(即血流动力学)将被增强。值得注意的是，该方法不需要用户指定阈值。函数hmrMotionCorrectCbsi用于执行此操作。

运动伪影减少算法(MARA)——MARA是一种校正方法，通过评估短时间窗口的振幅变化或标准差超过用户指定的阈值来检测伪影。发现包含伪影的时间窗被加宽，然后用平滑样条近似伪影，然后减去它以得到运动校正信号。函数hmrMotionArtifactByChannel用于检测运动时间段，而hmrMotionCorrectSpline用于校正。使用了Cooper等人(2012)的参数：tMotion=0.5，tMask=2，STDEVthresh=20，AMPthresh=0.5，pSpline=0.99。

目标主成分分析(tPCA)——第一个被广泛应用的去除fNIRS信号伪影的方法是PCA。该方法将fNIRS通道集提交到PCA（一种特征分解算法，将输入信号集转化为主成分的正交集）。然后排除那些占最大方差的成分，并从假定的非伪影成分重建信号。该方法的扩展版为目标主成分分析，是后来创建的，只在检测到伪影的短时间段上应用原始的PCA方法。这是使用hmrMotionCorrectPCArecurse函数实现的。使用的运动检测参数与MARA方法相同，原文使用到的其他参数包含nSV=0.97和maxIter=3。

峰度小波(kWavelet)——小波滤波是对输入信号进行离散小波变换，剔除幅值超过阈值的小波系数，再对信号进行重构的校正方法。该工作后来的扩展使用了小波系数的峰度来自动设置阈值。使用函数hmrMotionCorrectKurtosisWavelet，并将峰度阈值参数设置为原文推荐的3.3。

样条Savitzky-Golay(Spline-SG)——Spline-SG方法是一种多步骤校正过程，旨在通过减少用户提供的参数数量和增加一个单独的去除快速峰值的过程来扩展MARA。该方法首先使用Sobel滤波器提取信号的梯度，从中识别包含运动的时间段为超过1.5倍的四分位区间。然后，该方法尝试识别基线偏移和慢峰值，并使用MARA中的样条方法校正这些时间段。接着用Savitsky-Golay平滑滤波器去除快速尖峰。函数hmrMotionCorrectSplineSG使用的参数为：p=.99和FrameSize_sec=10。

激活。然后，按照如前所述(分析#1中的激活)的相同方式评估每种方法的激活性能。每次检验的AUC均方误差均小于10^-5。

实验验证

被试：23名典型发育的参与者参与了本次研究，年龄7-15岁(M=11.72岁，SD=2.13岁)。参与者排除标准包括：(1)用韦氏儿童智力量表(WISC-IV)或韦氏简化智力量表(WASI)测量的全量表智商低于80；(2)神经诊断(如癫痫)，以及基于父母的报告；(3)基于儿童和青少年症状量表-4R的精神病学诊断。

任务程序：实验环节包括一个7分钟的n-back任务。研究人员向参与者展示了一系列单一的辅音字母，并指示他们在呈现的当前字母与之前呈现的第n个字母相同时，用惯用手进行按键反应。n=0、1和2，分别是0-back、1-back和2-back block。0-back条件要求参与者只在看到“X”时进行按键反应，此时不会占用工作记忆。采用改进的拉丁方阵对刺激呈现的顺序进行伪随机化。每个block包含9个试次，每个试次持续3000ms，其中字母呈现500ms后，会有2500ms空屏。每个block为27s，block之间的间隔时间固定为14s，这样能够使血流动力学反应恢复到基线水平，然后给出5s的提示，即提示下一个block即将开始。在扫描之前，被试练习了每个任务条件的两个练习block。

成像程序：用双波长(690和830nm)连续波CW6(TechEn，Inc.，Milford，MA)成像系统记录光信号。以50Hz的采样率从探测器收集数据。56个光通道由16个光源和29个探测器组成，并安装在弹性纤维帽上(EASYCAP，Brain Products GmbH，Germany)。将帽子放置在头部模型上，连接EEG系统以确定标准10-20空间的光极坐标。然后使用AtlasViewer程序将光极和通道位置配准到Colin27大脑模板(图3)。

图3.探测器的布局。

分析

预处理。使用NIRS Brain AnalyzIR工具箱进行预处理和激活分析。然后将信号重采样至5Hz。

激活。通过创建工作记忆负荷依变回归来量化任务激活。如前所述，通过创建一个加权任务boxcar函数，其中每个条件(0-、1-和2-back)都有对应于其相对负荷(分别为1、2和3)的振幅。然后将加权的boxcar与正则血流动力学响应函数(HRF)进行卷积，并输入到一般线性模型。为了考虑血流动力学响应函数的变化，模型中包括了时间导数和离散导数。通过在模型中包含离散余弦变换项来解释低频漂移，最大频率为1/128Hz。通过Abdelnour和Huppert(2011)描述的基于正则化表面的图像重建，使用AtlasViewer在探测器配准过程中生成的灵敏度矩阵进行组水平分析。

结果

模拟#1：采样率和低通滤波的影响

模拟评估了采样率和低通滤波(LP)对TDDR性能的影响，发现当使用无预滤波的核心TDDR算法时，性能高度依赖于输入数据的采样率，以2Hz为最佳(图4)。然而，在所有高于1Hz的采样率下，添加低通预滤波步骤完全恢复了算法的性能。

图4.采样率和低通滤波对TDDR性能的影响。

模拟#2：方法比较

评估激活检测性能的ROC曲线如图5A所示。TDDR(0.775)的AUC均值最接近无运动伪影的数据(0.869)，大于CBSI(0.733)、MARA(0.563)、tPCA(0.591)、kWavelet(0.513)、spline-SG(0.652)和uncorrected(0.516)。各方法AUC值的均值和标准差如图5B所示。

图5.TDDR与其他校正方法的性能比较。

实验

每种方法对含氧血红蛋白数据产生的激活图谱如图6所示。在测试的方法中，与CBSI(3.02)、MARA(2.96)、tPCA(3.79)、kWavelet(3.96)、spline-SG(3.68)和未校正(3.67)相比，TDDR产生了最大的激活t统计量(4.88)。此外，与CBSI(903)、MARA(924)、tPCA(891)、kWavelet(935)、spline-SG(1153)和未校正(1560)相比，TDDR产生的正显著网格顶点数量最多(2399)(p＜.05)。

图6.不同方法激活结果的比较。

结论

本研究设计了TDDR运动校正方法。由于该方法基于鲁棒估计量，它没有调谐参数，可以有效地去除峰值和基线偏移。通过模拟和实证数据集的验证表明，该方法的性能优于目前的其他校正方法。这种方法增加了fNIRS应用于功能性神经成像的可行性，如应用于婴幼儿群体，以及那些因顺应性或交流限制而难以限制头部运动的智力障碍人群。

TDDR操作步骤（示例）

本研究提出的TDDR校正方法的操作步骤示例如下(该脚本由Robert Luke提供，基于python运行)。

import os
import mne

from mne.preprocessing.nirs import (optical_density,
                                    temporal_derivative_distribution_repair)

①加载数据

这里将使用fNIRS motor数据(该数据集来自麦考瑞大学)。重采样数据以使索引时间更快捷。然后，将数据转换为光密度，以对这些信号进行校正和绘制。

fnirs_data_folder = mne.datasets.fnirs_motor.data_path()
fnirs_cw_amplitude_dir = os.path.join(fnirs_data_folder, 'Participant-1')
raw_intensity = mne.io.read_raw_nirx(fnirs_cw_amplitude_dir, verbose=True)
raw_intensity.load_data().resample(3, npad="auto")
raw_od = optical_density(raw_intensity)
new_annotations = mne.Annotations([31, 187, 317], [8, 8, 8],
                                  ["Movement", "Movement", "Movement"])
raw_od.set_annotations(new_annotations)
raw_od.plot(n_channels=15, duration=400, show_scrollbars=False)

OUT：Loading /home/circleci/mne_data/MNE-fNIRS-motor-data/Participant-1

Reading 0 ... 23238 =      0.000 ...  2974.464 secs...

可以看到40、190和320s左右的运动中有一些小的伪影。但此数据相对干净，因此接下来将添加一些额外的伪影进行演示。

②向数据中添加伪影

NIRS数据中两种常见的伪影类型是尖峰和基线偏移。当一个人移动，并且光极相对于头皮移动然后返回其原始位置时，通常会出现尖峰。如果光极相对于头皮移动并且不返回其原始位置，则会发生基线偏移。这里，在数据中添加了100s的尖峰类型伪影和200s的基线偏移。

corrupted_data = raw_od.get_data()
corrupted_data[:, 298:302] = corrupted_data[:, 298:302] - 0.06
corrupted_data[:, 450:750] = corrupted_data[:, 450:750] + 0.03
corrupted_od = mne.io.RawArray(corrupted_data, raw_od.info,
                               first_samp=raw_od.first_samp)
new_annotations.append([95, 145, 245], [10, 10, 10],
                       ["Spike", "Baseline", "Baseline"])
corrupted_od.set_annotations(new_annotations)

corrupted_od.plot(n_channels=15, duration=400, show_scrollbars=False)

OUT：Creating RawArray with float64 data, n_channels=56, n_times=8924

Range : 0 ... 8923 =      0.000 ...  2974.333 secs

Ready.

③应用时间导数分布修复(TDDR)

这种方法无需任何用户提供的参数，即可校正基线偏移和尖峰伪影。

corrected_tddr = temporal_derivative_distribution_repair(corrupted_od)
corrected_tddr.plot(n_channels=15, duration=400, show_scrollbars=False)

我们可以从上面的数据中看到，虽然有一些残留的较小的伪影仍然存在，但原始数据中的伪影，以及引入的尖峰和偏移伪影在很大程度上被消除了。

注：需要TDDR工具包和示例脚本的小伙伴，可以在公众号后台回复“TDDR”获取。

参考文献：

Temporal Derivative Distribution Repair (TDDR): A motion correction method for fNIRS.

https://doi.org/10.1016/j.neuroimage.2018.09.025

小伙伴们点个“在看”，加(星标)关注茗创科技，将第一时间收到精彩内容推送哦～

更多精彩课程推荐

• 本文内容仅供参考，一切内容以英文原文为准。
  

• 非常欢迎大家留言、转载、收藏或分享~
• 本文来自微信公众号“茗创科技”。如需转载，请在“茗创科技”后台回复“转载”，并附上所需转载的文章标题以及您的ID。

MC_Brain

茗创科技工作室

觉得有帮助，欢迎转发收藏或者点个在看哦~

听说点在看的人SCI接收率都提升了18%呢！